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ABSTRACT 


Decoding the key dynamical processes that shape the Galactic disk structure is crucial for recon- 
structing the Milky Way’s evolution history. The second Gaia data release unveils a novel wave pattern 
in the Lz — (Vr) space, but its formation mechanism remains elusive due to the intricate nature of 
involved perturbations and the challenges in disentangling their effects. Utilizing the latest Gaia DR3 
data, we find that the Lz — (Vr) wave systematically shifts towards lower Lz for dynamically hotter 
stars. The amplitude of this phase shift between stars of different dynamical hotness (ALz) peaks 
around 2300kms~!kpc. To differentiate the role of different perturbations, we perform three sets of 
test particle simulations, wherein a satellite galaxy, corotating transient spiral arms, and a bar plus the 
corotating transient spiral arms act as the sole perturber, respectively. Under the satellite impact, the 
phase shift amplitude decreases towards higher Lz, which we interpret through a toy model of radial 
phase mixing. While the corotating transient spiral arms do not generate an azimuthally universal 
phase shift variation pattern, combining the bar and spirals generates a characteristic Alz peak at 2:1 
Outer Lindblad Resonance, qualitatively resembling the observed feature. Therefore, the Lz — (Vp) is 
more likely of internal origin. Furthermore, linking the ALz peak to the 2:1 Lindblad resonance offers 
a novel approach to estimating the pattern speed of the Galactic Bar, supporting the long/slow bar 
model. 
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1. INTRODUCTION 


The Accurate astrometric information provided by the 
Gaia mission has revolutionized the field of Milky Way 
dynamics. Three-dimensional positions and velocities 
of millions of stars provided by synergies between Gaia 
DR2 (Gaia Collaboration et al. 2018) and large spec- 
troscopic surveys have led to the discovery of a series of 
phase space substructures in the Galactic disk, including 
diagonal ridges in the R— Va plane (Antoja et al. 2018; 
Kawata et al. 2018; Ramos et al. 2018). The poten- 
tial contributing perturbers for generating these R — Vz 
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ridges and associated velocity substructures include the 
Galactic bar (Mühlbauer & Dehnen 2003; Chakrabarty 
2007; Antoja et al. 2018; Monari et al. 2019; Fragkoudi 
et al. 2019), spiral arms (Quillen et al. 2018; Hunt et al. 
2018, 2019; Khanna et al. 2019), and the Sagittarius-like 
satellite (Minchev et al. 2009; Gómez et al. 2012; La- 
porte et al. 2019; Khanna et al. 2019). Understanding 
the physical origin of these ridges is crucial for recon- 
structing the Milky Way’s dynamical evolution history. 

However, a consensus on the formation mechanisms 
of these ridge-like structures remains elusive, owing to 
the complexities associated with various perturbations 
and the challenges involved in disentangling their ef- 
fects. Test particle simulations conducted by Hunt et al. 
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(2019) have demonstrated that the combination of a 
bar with arbitrary pattern speeds and transient spirals, 
can qualitatively reproduce the observed R — Vg ridges, 
thereby making it exceedingly difficult to isolate the in- 
dividual effects of these perturbers. According to the 
N-body simulation sets of Khanna et al. (2019), both 
transient spiral arms excited in isolation, and a satellite 
perturber at the mass scale of 101° Mọ generate ridges 
qualitatively similar to the observations in the density 
and (Vz) map. This complexity necessitates the identi- 
fication of new discriminating features. 

Ridge-like structures exhibit a strong connection with 
radial motion, as evidenced by the mean radial velocity 
(Vr) map in the R — Ma plane (Fragkoudi et al. 2019; 
Hunt et al. 2019; Khanna et al. 2019; Wang et al. 2020). 
Notably, these structures align approximately along the 
lines of constant angular momentum (Lz = Rx V4), mir- 
roring the morphology of ridges in the number density 
map. Utilizing orbit integration in an N body poten- 
tial, Fragkoudi et al. (2019) demonstrated that a long 
R — V; ridge, accompanied by a change in the (Vp) 
direction, could pinpoint the location of the 2:1 outer 
Lindblad resonance, reinforcing the physical connection 
between ridges in these two projections. Consequently, 
the (Vp) corrugation binned in angular momentum can 
be viewed as the one-dimensional projection of the 2D 
ridges, simplifying comparative analysis while retaining 
the essential physical information. 

Upon first discovery by Friske & Schönrich (2019), 
this Lz — (Vr) wave displays systematic displacement 
towards lower Lz for stars with higher vertical energy 
(Ez), suggesting a phase shift among the wave pat- 
tern of stars with different dynamical hotness. Friske & 
Schénrich (2019) attributed this feature to orbital reso- 
nances, which raises an interesting question on the exis- 
tence of such phase shift when subject to other pertur- 
bations like the satellite pericenter passage. Moreover, 
it remains unknown whether the variation in the phase 
shift amplitude with Lz encodes information about its 
origin. 

Previous works have shown that phase space substruc- 
tures could vary with dynamical hotness. In the (Vp) 
color-coded R — Vy space, Wang et al. (2020) found 
some ridges varying among populations of different stel- 
lar ages whereas others do not. Since stellar age is a 
proxy of dynamical hotness, interpreting this dichotomy 
from the dynamical perspective may help uncover its ori- 
gin. In the vertical direction, Li & Shen (2020) found 
the Z — Vz phase spiral becomes less prominent or even 
absent for those stars with higher Jz (dynamically hot- 
ter in the radial direction). Analogously, an analytic 
model of Laporte et al. (2020) also demonstrated that 


the R— Vy ridges generated by bar resonances are more 
prominent for the dynamically colder population. 

To explore the connection between ridges and spiral 
arms arising from a single satellite impact, Antoja et al. 
(2022) developed an analytical model and found that the 
R — Vz ridges exhibit a V-shaped morphology in both 
test particle and N body simulations. However, an intu- 
itive explanation for the formation of such morphology 
remains lacking. Furthermore, an undulating Lz — (Vr) 
wave emerges simultaneously, with its frequency increas- 
ing during the phase mixing. Fourier transform of the 
wave reveals two frequency peaks, attributed to pertur- 
bations occurring less than 0.4 Gyr ago and 0.7-1.8 Gyr 
ago, respectively. However, their conclusions only hold 
if the satellite is the sole perturber. Both Antoja et al. 
(2018) and Khanna et al. (2019) used toy models of 
winding spirals to mimic the R — Vy ridges and found 
qualitative agreement with the observed R — Vy den- 
sity map, without accounting for its correlation with the 
(Vr) map. A deeper understanding of the role of other 
perturbing mechanisms in generating the Lz— (Vpr) wave 
is still lacking. 

We present a novel perspective to ascertain the origin 
of the Lz—(Vp) wave by analyzing its dependence on dy- 
namical hotness. We quantify this dependence through 
the phase shift between waves of varying dynamical hot- 
ness. In Section 2, we propose a toy model of radial 
phase mixing to intuitively understand the formation of 
the Lz — (Vr) wave after the external satellite pertur- 
bation and use it to make a theoretical prediction on 
phase shift variation trend. Observational analysis of 
the phase shift variation pattern is in Section 3.2. We 
conduct two groups of test particle simulations in Sec- 
tion 3 to differentiate the role of internal and external 
perturbers in generating the phase shift variation pat- 
tern. After presenting the simulation results in Sections 
4.1 and 4.2, we discuss the caveats and future directions 
in Section 5, and summarize our main findings in Section 


6. 


2. TOY MODEL OF RADIAL PHASE MIXING 


In this section, we introduce a simple toy model 
of radial phase mixing to elucidate the formation of 
phase space substructures following external perturba- 
tion. Additionally, we employ this model to qualita- 
tively predict the dependence of the Lz — (Vr) wave 
on dynamical hotness. We represent each Lz bin with 
a single particle orbit, and the phase space coordinates 
of different orbits (corresponding to different Lz bins) 
are stitched together to depict regions of relatively high 
phase space density. The radial velocity Vz of the orbit 
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Figure 1. Illustration of the toy model of radial phase mixing presented in Section 2. The left panel maps the distribution 
of particles in the Qpr — Or and the R — Va space, both color-coded by (VR). The right panel illustrates the Lz — (Vr) wave 
signal of the dynamically cold (Jz < 8kms™* kpc, blue) and hot (Jz > 23kms~' kpc, red) populations at the time of 250 and 
400 Myr. As phase mixing proceeds, the characteristic wavelength of the Lz — (Vr) decreases. The wave of the dynamically 
hot population displays a systematic phase shift towards lower Lz compared to the cold one. The phase shift amplitude (ALz, 
blue points) increases with decreasing Lz, as predicted by our toy model. 


represents the mean radial velocity (Vr) of its corre- 
sponding Lz range. 

To give the toy model predictive power in a realistic 
sense, we utilize the same initial condition as the test 
particle simulations in Section 3.3 and define the dy- 
namically cold and hot sub-sample in the same manner. 
We also calculate the mean values of radial frequency 
Qpr for each Lz bin using agama (Vasiliev 2019). The 
OR value is given by the equation: 


Orp=ORXT+4R 0 (1) 
where T is the perturbation time. For simplicity, we 
assume both subpopulations follow the same distribu- 
tion in the Qr — OR plane. While this assumption is not 
strictly correct, it has a negligible impact on the key 
features of our prediction. By requiring all particles to 


start at the same radial phase (0r, = 0), we implicitly 
assume that all stars receive the impulsive perturbation 
from a distant perturber. Then we utilize epicycle ap- 
proximation to compute the phase space coordinates for 
the particles using the following equations: 
Vr = X cos r, R = Rg + X sindR/OQR (2) 
(See Binney & Tremaine 2008, for detailed derivation). 
Upon adopting this equation, we also implicitly assume 
the particle orbits remain regular after the satellite per- 
turbation. Then we map their distribution in phase 
space projections like the Lz — (Vpr) and the R — Vy 
to understand the structural correlation between them. 
The numerical prescription is as follows. We 
choose particles with angular momentum Lz € 
[800, 3000] kms~!kpc. The corresponding guiding ra- 
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dius Ry is inferred from the rotation curve of the adopted 
Milky Way potential model MWPotential2014 (Bovy 
2015). The radial oscillation amplitude X is in the form 
of X = Ax Ry. Changing the form of the dependence 
of X on guiding radius R, has a negligible effect on the 
main results. The value of A is set at 1.2 to ensure the 
(Vr) amplitude is at the same order of magnitude as 
both the observations and simulation results. 

As the upper left panel of Figure 1 shows, a series of 
parallel lines, each with the same slope set by perturba- 
tion time, emerges in the Qpr — OR plane. In the bottom 
left panel, ridge lines color-coded by Vpr form connected 
V-shaped streaks. The adjacent branches of the “V” 
shape display opposite signs of (Vr). Ridge lines reach 
a turning point when (Mo) changes sign and Oo reaches 
0 (27) or m. These two features are qualitatively consis- 
tent with the N-body simulation of the interaction be- 
tween an Sgr-like satellite and a Milky Way-like galaxy, 
as shown in Fig.18 of Antoja et al. (2022). 

The Lz — (Vr) wave also emerges, where (Vr) = 0 
corresponds to a turning point of a V-shaped streak in 
the R— Vy plane. The spacing between the adjacent 
zero radial velocity points becomes wider at higher Lz 
because of the shallower slope of the Lz — Qr curve. 
The phase shift between waves of different dynamical 
hotness arises from differences in mean Ur. Since we 
assume both populations follow the same (Ur — Op line 
series, this frequency difference leads to the phase dif- 
ference in Or. This phase difference increases towards 
the lower Lz range, where the difference in dynamical 
hotness (or mean Q,)increases. As shown in the right 
panel of Figure 1, the phase difference in Or increases, 
resulting in a shift in extrema location, with its ampli- 
tude (ALz, which will be defined in the next section) 
increasing with decreasing Lz. As time increases, the 
slopes of the QR — Op line series become steeper, which 
causes the characteristic wavelength (or frequency) of 
the Lz — (Vr) wave to decrease (or increase), as the 
right column of Figure 1 demonstrates. Antoja et al. 
(2022) utilized this property to date the perturbation. 
From Equation 1, the phase difference in Or increases in 
the meantime. However, the characteristic wavelength 
of the Lz — (VR) wave also decreases, which cancels out 
the effect of increasing phase difference in Or and main- 
tains a roughly constant phase shift amplitude. Since 
this phase shift accumulation with Dz is due to the dif- 
ference in Qpr, its maximal value occurs in the low Lz 
end, rather than at the high Lz end where the external 
satellite perturbation is strongest. 


3. METHODOLOGY 
3.1. Phase Shift Measurement 


We define the phase shift amplitude as the differ- 
ence between the locations of local extrema, ALz = 
LZcold — Lz,not- The subscripts correspond to the dy- 
namically cold and hot sub-samples of the whole distri- 
bution, which is dissected based on the values of Jz. To 
analyze the variation trend of ALz with Lz (referred 
to as Lz variation trend for brevity in the following 
text), we define the mean Lz location of the extrema 
as Lz = (Lz coia + Lznot)/2. Strictly speaking, quan- 
tifying the actual phase shift requires dividing ALz by 
the characteristic wavelength of the wave. Nevertheless, 
we adopt the value of ALz as a proxy for phase shift be- 
cause the ever-changing shape of the wave across differ- 
ent Lz ranges severely complicates the task of extract- 
ing characteristic wavelengths and may introduce addi- 
tional systematic bias. We focus on the Lz — (Vp) wave 
pattern in the range of Lz € [600,3000] kms~! kpc, di- 
vided into 165 equally-spaced bins. First, we smooth 
the curve with a Gaussian kernel with a size of (o =) 4 
Lz bins (58.2 kms~'kpc) to mitigate the effect of 
small-scale noise. Visual inspection of the smoothed 
wave signal ensures that this process does not gener- 
ate pseudo-oscillations that may compromise the phase 
shift measurement. Then, we adjust the parameters 
of scipy.signal.find_peaks function (i.e. distance, 
width etc.) to find the suitable parameter set capable of 
identifying all noticeable extrema points of the Lz— (Vp) 
wave pattern, and apply the same settings to all analy- 
ses with simulations. We discard those extrema points 
that have no counterparts in the other population. Fur- 
thermore, we set the distance parameter at 11 to avoid 
adjacent extrema points being too close to each other. 
The upper bound of the width parameter is also set 
at a value high enough to identify extrema points at 
higher Lz where the wavelength becomes longer in that 
range in the case of external satellite perturbation. We 
only include an extrema point when its prominence pa- 
rameter is greater than 1 to mitigate the effect of small 
amplitude oscillation. 

With the above setup, we divide each simulation snap- 
shot into eight equally spaced azimuthal ranges and ex- 
tract extrema points from the “cold” and “hot” waves. 
We compare the Lz variation trends in different az- 
imuths to conclude a universal variation pattern. If 
there is none, we try to unveil the ALz variation feature 
existing in most azimuthal ranges, which may also be 
valuable for discriminative purposes. Due to the subtle 
and discrete nature of the measurable extrema points, 
conclusions drawn from the analysis on the ALz vari- 
ation trend are only reliable in the qualitative sense, 
and any quantitative conclusion should be treated with 
great caution. A quantitative match of the ALz vari- 


ation curve between simulation results and observation 
data is beyond our capability given the simplicity of our 
simulation configuration. 


3.2. Data 


Among 33 million stars with line-of-sight velocity 
measurement from Gaia DR3, we obtain 26,611,026 
sources meeting the criteria for reliable StarHorse (An- 
ders et al. 2022) distances, with fidelity > 0.5 and 
sh-outflag==” 0000”. StarHorse is a Bayesian code 
that leverages astrometric, photometric, and spectro- 
scopic data from multiple surveys to derive the cu- 
mulative distribution function of astrophysical param- 
eters, including the Heliocentric distance (whose me- 
dian is dist50 column). We adopt a distance of 
8.275 kpc between the Sun and the Galactic center 
(GRAVITY Collaboration et al. 2021), with the Sun 
situated 20.8 pc above the Galactic midplane (Ben- 
nett & Bovy 2019). For the motion of Sgr A*, we 
use a radial velocity of —8.4kms~! (GRAVITY Col- 
laboration et al. 2021) and a proper motion in the 
ICRS frame picrs = (—3.16, —5.59) masyr~ (Reid & 
Brunthaler 2020). Combining these measurements 
yields the adopted solar peculiar velocity components 
(Uo, Vo, Wo) = (8.7, 251.5, 8.4) kms~!. Our cut in the 
azimuthal range |¢ — ¢o| < 0.2rad gives a sample size 
of 19,279,240. Using the “Stackel Fudge” method (Bin- 
ney 2012) incorporated in the agama package (Vasiliev 
2019), we calculate the action-angle-frequency quanti- 
ties for the entire sample, employing MWPotential2014 
(Bovy 2015) as the Milky Way potential model. We 
choose not to consider the selection function effect in 
our study due to its minor impact on the mean velocity 
map. 

As depicted by the solid green line of Figure 2, the 
overall shape of the Lz — (Vr) wave for the whole sam- 
ple is consistent with previous works(Friske & Schénrich 
2019; Antoja et al. 2022). To examine the morphologi- 
cal variation of the Lz — (Vp) wave with dynamical hot- 
ness, we categorize stars based on their vertical action 
Jz: stars with Jz < 3kms~! kpc are considered dynam- 
ically “cold,” while those with Jz > 12kms~!kpc are 
considered dynamically “hot”. This division guarantees 
a sufficient difference in dynamical behavior between the 
two subsamples, allowing for measurable phase shifts in 
their wave patterns, while also ensuring that the mor- 
phological difference between the wave patterns is not 
substantial enough to introduce systematic bias in the 
phase shift measurement. The investigated Lz range is 
narrower ([800, 2800] kms~! kpc), beyond which the ob- 
servational uncertainty is too great for robust phase shift 
measurement. The parameters of the peak-finding pro- 
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Figure 2. Shape and phase shift between particles with dif- 
ferent orbit hotness of the Lz — (Vr) wave in the Gaia DR3 
data. The green line is the Lz — (Vr) wave of the whole 
sample. The Lz — (Vr) wave pattern composed of dynami- 
cally cold (Jz < 3kms™' kpc) and hot (Jz > 12kms~' kpc) 
are depicted in red and blue lines, respectively. The brown 
dashed line displays the phase shift amplitude (ALz) varia- 
tion pattern which peaks at around 2300 kms~! kpc. 


cedure are the same as described in Section 3.1 except 
for the Gaussian kernel size (o = 2). The extrema points 
of the “hot” wave systematically shift towards lower Lz 
compared to the “cold” wave, qualitatively consistent 
with the results of Friske & Schonrich (2019) who used 
vertical energy Ez as a proxy for dynamical hotness. 
Both the “hot” and “cold” waves exhibit similar shapes, 
except at the high Lz end near 2800kms~!kpc. Due to 
the larger measurement uncertainty of the data, whether 
the two waves are truly phase-aligned at the high Lz end 
is debatable. As illustrated by the brown dashed line of 
Figure 2, the phase shift amplitude (denoted as ALz) 
is largest at the extrema point near 2300 kms~!kpc, 
which presents a prominent peak in the ALz variation 
trend. 

The phase shift ALz variation trend in observations 
raises some interesting questions: does this trend differ 
among different types of perturbations? If they do, can 
they be used to provide constraints on the dynamical 
evolution history of the Milky Way disk? Can we find 
a way to qualitatively understand the physical origin of 
this trend? 


3.3. Test particle simulations 


We now turn to test particle simulations involving in- 
ternal and external perturbations, to examine the dif- 
ference in the resulting phase shift variation pattern of 
the Lz — (Vr) wave. Compared to N-body simulations, 
test particle simulations neglect self-gravity, but do al- 
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low us to run more particles and test the effect of varying 
specific parameters on the orbit evolution. 

We sample eight million disk particles from the 
quasiisothermal df (Carlberg & Sellwood 1985; Bin- 
ney 2010; Binney & McMillan 2011) implemented 
in agama (Vasiliev 2019). These particles are dis- 
tributed using a radial scale length of 3 kpc within 
the Milky Way potential model, MWPotential2014 
(Bovy 2015). The central velocity dispersion 
(or|R=0, 0z|R=0) = (90,110) km es" gives a disk that is 
dynamically colder than the actual Milky Way disk. The 
scale lengths of the og and oz profiles are 6 and 7 kpc. 
The Lz—(Vr) wave signal is more prominent, which sim- 
plifies the task to study the Alz variation trend with 
Lz. To investigate the impact of dynamical hotness on 
the Lz—(Vp) wave morphology, we define those particles 
with Jz < 8kms~!kpc as the dynamically cold” popu- 
lation, and those with Jz > 23kms~! kpc as the dynam- 
ically "hot" population. Under such division, the par- 
ticle number of the cold and hot populations is roughly 
the same. 

The external perturber (Milky Way satellite) is a 
2.5 x 101° Mo Plummer sphere with a half-mass ra- 
dius of 3 kpc. Its position 7 = (4,8,18)kpe 
and velocity ¢ = (—339,—44,76)kms~! at the first 
pericenter passage are from the El model of de 
la Vega et al. (2015). Backward orbit integration 
for 1 Gyr in the MWPotential2014 potential using 
galpy’s ChandrasekharDynamicalFrictionForce rou- 
tine to account for dynamical friction provides the initial 
condition for our simulation. As Figure 3 shows, the per- 
turber only impacts the Galactic disk once during the 
first 2 Gyr of the simulation. At 2 Gyr, we reduce the 
mass and radius of the satellite to 1 x 101° Mo and 1 kpe 
to mimic the mass loss after the first pericenter passage. 
The second pericenter passage occurs at 3.1 Gyr when 
the satellite crosses the disk plane at R = 15kpc with 
Vz % 300kms~!. The total integration time is set to 4 
Gyr to cover these two pericenter passages. 

Using test particle simulations, we also investigate 
the influence of internal perturbations with the com- 
bination of a steadily rotating bar and two-armed tran- 
sient spirals. The spiral arm potential reaches its maxi- 
mal amplitude at 1200 Myr. To differentiate the role 
of different resonances in generating the ALz varia- 
tion trend, we complement this simulation with a test 
where the transient spirals act as the sole perturber. 
In this test, the spiral potential has a higher ampli- 
tude and reaches the maximal value at the start. The 
2-armed transient spirals follow the default model de- 
scribed in Section 2.2 of Hunt et al. (2018) which fol- 
lows the potential form given by Cox & Gómez (2002). 
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Figure 3. The orbit of the Sgr-like satellite in our test 
particle simulations, following the prescription in Section 3.3. 
The circles and crosses represent the orbits of the satellite 
before and after reducing its mass from 2.5 x 10'°Mo to 1x 
10'°Mo. Each point is color-coded with time at an interval 
of 50 Myr. The bold gray line represents the Milky Way disk. 


The winding and decay of the spiral arms are con- 
figured by the CorotatingRotationWrapperPotential 
and GaussianAmplitudeWrapperPotential routine of 
galpy. The morphological parameters of the spiral 
model are N = 2, Rs = 0.3, Cn = 1, H = 0.125, fsp = 12° 
(sp represents the pitch angle). The bar is modeled 
as the 3D generalization of the Dehnen bar potential 
(Dehnen 2000; Monari et al. 2016) with the same pat- 
tern speed (40 kms~!kpc~', invariable with time) and 
bar length (4.5 kpc) as the “Fiducial” model presented 
in Trick et al. (2021). The integration time is 2 Gyr for 
the case with a bar and 0.4 Gyr without a bar. Par- 
ticle orbits are integrated using the galpy code (Bovy 
2015). Previous literature shows that a steadily rotat- 
ing bar can only produce at most two prominent R — Vg 
ridges (Antoja et al. 2018; Hunt et al. 2018), we decide 
not to show the same results of the bar here to avoid 
redundancy. Generating more ridges requires higher or- 
der Fourier modes from the bar potential (Monari et al. 
2019), which is beyond the scope of this work. In addi- 
tion, the number of measurable extrema points would be 
too limited to conclude on a clear ALz variation trend. 
We exclude stars with an initial radius smaller than 1.5 
kpc to save integration time. This will have little impact 
on the simulation results since the number of discarded 
particles capable of entering the investigated Dz range 
is negligible. 


4. SIMULATION RESULTS 
4.1. External Satellite Perturbation 
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Figure 4. Temporal evolution of the modeled Galactic disk after the external satellite perturbation. The upper and lower rows 
are number density and the mean radial velocity (Vg) map in the X — Y plane. The evolution times of the four snapshots are 
1400, 1600, 2000, and 3400 Myr from left to right. The first and second pericenter passages are at around 1 and 3.1 Gyr. As 
phase mixing proceeds, the spiral arms become more tightly wound. 


Figure 4 illustrates the evolution of the MW-like disk 
after the satellite pericenter passage. Tidally induced 
spiral wraps are visible in the X — Y projection of sur- 
face density and mean radial velocity (Vr). As time pro- 
gresses, the spirals become more tightly wound. (Vp) 
of the innermost region is close to zero due to swifter 
phase mixing which expands in spatial extent with 
time. Simultaneously, the characteristic wavelength of 
the Lz — (Vr) wave shown in Figure 5 increases with 
Lz, consistent with our toy model in Section 2 and the 
model of Antoja et al. (2022). After the second pericen- 
ter passage, the spiral pattern exhibits more complex 
morphology, as shown in the rightmost column of Fig- 
ure 4. We do not show snapshots taken at the earlier 
time after the satellite’s pericenter passage because the 
number of measurable extrema points is too few to re- 
veal a clear trend. 

Two trends are visible in the ALz — Lz plot of Fig- 
ure 5. At a fixed time, Alz decreases as Lz(Qpr) in- 
creases(decreases) but rarely becomes negative in the 
lower Lz range, which is consistent with the prediction 
of our toy model presented in Figure 1. The amplitude 
of ALz is also at the same order of magnitude as our 
prediction. The relatively large ALz at the high Lz 


end may be related to the irregular morphology of the 
orbits strongly perturbed by the satellite, which inval- 
idates the usage of epicycle approximation in our toy 
model. This suggests that our toy model does not apply 
to the outermost part of the disc. Besides that, at a 
fixed Lz range, the phase shift amplitude ALz displays 
negligible variation with time. This also aligns with our 
theoretical prediction in Section 2. Our results unveil a 
novel formation channel for the phase shift in addition 
to the orbital resonance suggested by Friske & Schénrich 
(2019). 

As the lower right panel of Figure 5 demonstrates, 
the Lz — (Vr) wave signal exhibits less resemblance to 
the sinusoidal wave after the second pericenter passage. 
However, the decreasing A Lz trend persists in the lower 
Lz range. The temporal variation of A Lz is also negligi- 
ble. Therefore, the same ALz variation trends described 
above persist as long as the external satellite remains 
the dominant perturber, regardless of the wave shape or 
the number of pericenter passages. Despite our simple 
phase measurement, the clarity of these trends validates 
the physical process described in Section 2. 

The minor variation of phase shift with time has im- 
portant implications on two fronts. Firstly, it renders 
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Figure 5. Phase shift amplitude (ALz) variation trend versus Lz after the satellite pericenter passage. The plots in the upper 
and lower panel come from the simulation data at 1400, 1600, 2000, and 3400 Myr. The dashed lines in the upper row illustrate 
ALz — Lz plot in 8 different azimuthal ranges. In the lower row, the red and blue solid lines display the Lz — (Vr) wave 
signal for the dynamically hot and (Jz > 23kms~*kpc) cold (Jz < 8kms~‘*kpc) subpopulations. Within the lower Lz range, 
we see a monotonic decreasing trend of ALz with increasing Lz in most azimuthal ranges for all snapshots, consistent with the 


prediction of our toy model developed in Section 2. 


the utilization of ALz variation pattern to date pertur- 
bation infeasible. Fourier analysis technique proposed 
by Antoja et al. (2022) is more suitable for this task. 
On the other hand, this constancy suggests that this 
unique feature can aid in ascertaining the wave’s origin 
regardless of the evolutionary stage. 


4.2. Internal Perturbers: Bar & Transient Spiral 


In the case of the corotating transient spirals, the 
phase shift amplitude AL exhibits no universal trends 
with Lz among different azimuthal ranges. As Figure 6 
demonstrates, within the majority of azimuthal ranges, 
the phase shift amplitude ALz displays irregular oscilla- 
tion, with an amplitude between 30 and 80 kms~! kpc. 
We are not able to observe and ALz variation feature 
shared by both snapshots. Interestingly, the decreasing 
trend with Lz resembling the trend produced by exter- 
nal satellite perturbation exists for very few azimuthal 
ranges of the two snapshots. In such cases, the slope of 
the ALz curve is shallower than those generated from 
the satellite impact, which still allows us to distinguish 
from the ALz variation trend of external perturbation. 

When combining the perturbations of the bar and 
transient spirals, the disk evolution is shown in Figure 7. 
The central depression in density is the result of exclud- 
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Figure 6. The temporal evolution of the modeled Galactic 
Disk under the perturbation of 2-armed corotating transient 
spirals in the same manner as Figure 5. The time is at 200 
and 300 Myr. 


ing the innermost particles from orbit integration. Be- 
fore the transient spiral starts growing (T< 1000 Myr), 
the butterfly pattern in (Ve), characteristic of a bar is 
visible. The signs of (Vr) change at Lindblad resonances 
(Binney & Tremaine 2008). When the spiral potential 
reaches maximal amplitude at 1200 Myr, the two-armed 
spiral pattern becomes most prominent in both the den- 
sity and (Vp) maps. The amplitude of (Vp) decreases as 
the spirals wrap up and decay. The bar resumes domi- 
nation as the time approaches the end of our simulation. 

As Figure 8 illustrates, after the spirals become ac- 
tive, although the overall ALz pattern varies drastically 
with time, a common feature in the phase shift variation 
pattern emerges across most of the azimuthal ranges: 
a characteristic peak with timely variable prominence. 
Interestingly, the location of ALz peak coincides with 
the location of 2:1 outer Lindblad resonance (2:1 OLR, 
marked by a vertical blue dashed line) in most azimuthal 
ranges. In a minority of the azimuthal ranges, the Al: 
peak is less prominent or absent, which makes the AL z 
variation trend qualitatively similar to those produced 
by the transient spirals. Despite the caveats mentioned 
above, the sufficiently prominent ALz peak could pro- 
vide the most likely range of 2:1 OLR. In a case where 
we would not know the true value of pattern speed (in, 
we could use the Lz position of the 2:1 OLR peak and 
use it to infer the Q, (while assuming a rotation curve). 

Utilizing test particle simulations with bar as the sole 
perturber, Trick et al. (2021) (as also previous litera- 
ture like Mühlbauer & Dehnen 2003) proposed that the 
(Vr) map in Lz — Jp action space displaying change in 
sign is the signature of 2:1 OLR. However, owing to the 
interference of transient spiral perturbations, the sign 
change in the Lz — (Vp) wave pattern also occurs in Lz 
ranges other than specific types of resonances, as demon- 
strated by Hunt et al. (2019) in the R — Ma projection. 
Among the four candidates of 2:1 OLR, the “Hat” and 
“Sirius” moving groups are the two closest to the ALz 
peak. This gives a pattern speed ranging between 33 and 
41 kms~!kpe~!, which favors the long/slow bar model. 
The above constraint is consistent with previous work 
using the solar neighborhood kinematics to indirectly 
measure the bar pattern speed (Chiba & Schönrich 2021; 
Trick 2022). 

Due to the fundamental difference between internal 
and external perturbations, our toy model of radial 
phase mixing is unable to account for the ALz variation 
trends generated by the bar and spirals. One perspective 
that may help account is the dependence of resonance 
lines on dynamical hotness (Jz). As illustrated in Figure 
10 of Trick et al. (2021), the location of bar resonance 
(so-called ” ARL”, axisymmetric resonance line) shifts 


9 


towards lower Lz for stars of higher Jz. In the region 
of the bar’s Outer Lindblad resonance, the majority of 
stars are forced towards higher Lz, and the changes of 
Lz (dLz[Lz,0]) are smaller than corotation resonance. 
Under corotation resonance, radial(Lz) migration can 
proceed in both directions, which could partially cancel 
out the effect of ARL shift, which is not the case for 
Lindblad resonances. Therefore, the bar’s outer Lind- 
blad resonance can generate a ALz larger than the coro- 
tating transient spirals and leave a distinct ALz peak in 
the presence of the spirals. No prominent peaks emerge 
under corotating transient spiral perturbation since the 
orbit dynamics are dictated by corotation resonance in 
the whole Lz range. 

To briefly summarize, the Lz — (Vr) wave pattern 
displays different dependence on dynamical hotness de- 
pending on the nature of dominant perturbations. Uti- 
lizing the amplitude of phase shift (AL zl between the 
wave patterns of different dynamical hotness, we could 
quantify this dependence and associate its variation 
trend versus Lz with particular perturbation types. The 
external satellite generates a qualitative trend of de- 
creasing ALz in the lower Lz range in almost all az- 
imuths during the whole time. Such universality does 
not exist in the case of corotating transient spirals. 
Under the joint perturbation of the bar and the tran- 
sient spiral, ALz displays a characteristic peak around 
2:1 outer Lindblad resonance, which qualitatively re- 
produces the observed ALz feature. Controlled tests 
demonstrate that the combination of bar and corotat- 
ing transient spirals are the most likely perturbations 
among the three investigated cases. 


5. DISCUSSION 
5.1. Limitations 


Despite capturing the essential physical processes of 
radial phase mixing, our toy model presented in Section 
2 is not exempt from limitations. Firstly, it outlines 
the backbone of the phase space substructure but falls 
short of providing a detailed phase space density distri- 
bution. As the system becomes fully phase-mixed in the 
lower Lz range, the (Vr) amplitude should be close to 
zero, which is not demonstrated in the toy model. Ne- 
glecting the azimuthal dependence of the external per- 
turbation prohibits our toy model from describing the 
azimuthal phase mixing process as the model of Antoja 
et al. (2022) does. Employing the epicycle approxima- 
tion, our model can only account for the (Ur difference 
between cold and hot orbits, prohibiting us from ex- 
plaining the larger ALz at the high Lz end. 

It is also imperative to be cautious while interpret- 
ing the ALz variation trends generated by test parti- 
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Figure 7. The temporal evolution of the modeled Galactic Disk under the joint perturbation of corotating transient spiral 
and steadily rotating bar demonstrated in the same manner as Figure 4. The transient spiral potential reaches its maximal 
amplitude at 1200 Myr. 
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Figure 8. Phase shift (ALz) variation trend is shown in the same manner as Figure 5. Four snapshots are from simulation 
data at 800, 1200, 1500, and 1800 Myr. The vertical dashed lines mark the locations of corotation resonance (brown) and 2:1 
outer Lindblad resonance (blue). 


cle simulations, given the absence of self-gravity effects. 
Future investigations employing tailored N-body simu- 
lations are required to validate our main results. 

We should note that the ALz variation trend can only 
differentiate different perturbations in the probabilistic 
sense because not all azimuthal ranges exhibit the com- 
mon AL 7 variation features unique to specific perturba- 
tions. Under satellite impact, the decreasing AL z trend 
predicted by our toy model is less pronounced in some 
of the azimuthal ranges. In the case of bar plus spiral 
arms, there is no guarantee that the 2:1 OLR generates 
the characteristic ALz peak in all azimuthal ranges, ei- 
ther. Therefore, broader azimuthal coverage enabled by 
future surveys is essential for robust diagnostics. 


5.2. Comparisons with Previous Works 


Friske & Schénrich (2019) measured the phase shift by 
fitting the observed data with the sinusoidal function 
and comparing the best-fit phase angle between each 
other. The unit of their phase shift is the degree, not 
the angular momentum used (kms~! kpc) in our work. 
While their approach is advantageous in utilizing com- 
plete wave information, it is susceptible to systematic 
errors due to its model-dependent nature. During fit- 
ting, they set the wavelength of the sine function to fixed 
values, essentially making their approach akin to peak 
finding in principle. Our method does not assume any 
waveform, but the results of peak finding could vary with 
the width of Lz bins. Both approaches have the issue of 
defining the phase shift when the wave shapes deviate 
significantly from the sinusoidal shape. New mathemat- 
ical tools that can robustly measure the phase shift in 
such non-sine-like waves are indispensable for future re- 
search. 


5.3. Implications For Gaia DR Observation 


The comparison between observational data and our 
test particle simulations suggests that the Lz — (Vr) 
wave is more likely to originate from internal perturba- 
tions rather than external satellites. In the case of exter- 
nal perturbation, extrema points with the most promi- 
nent ALz are located at both ends of the investigated 
Lz range, which disagrees with the key observational 
feature. In contrast, the characteristic ALz peak re- 
lated to the bar’s 2:1 outer Lindblad resonance persists 
regardless of the evolution phase of the corotating tran- 
sient spiral arms. Therefore, if the observed phase shift 
variation trend does represent the trends within most of 
the azimuthal ranges, the combination of the bar and 
the transient spiral is more likely to be the dominant 
driving force. It is important to note that this conclu- 
sion is not definitive, as other scenarios like a deceler- 
ating bar (Chiba et al. 2021) have not been explored in 
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this work. Nevertheless, our study provides a novel ap- 
proach to differentiating the roles of internal and exter- 
nal perturbations in sculpting the observed phase space 
substructures. 

If we accept the bar and the transient spiral arms as 
the dominant perturbers, the peak of the phase shift 
amplitude AL is associated with 2:1 OLR. This ALz 
diagnostic could pave a new way to constrain the pattern 
speed of the Galactic bar. However, the uncertainty 
is greater than conventional methods because we can 
not exactly pinpoint the resonance location from several 
extrema points. Nevertheless, it provides a viable way 
to detect the bar resonance signature in the presence 
of spirals, which could help break the degeneracy in the 
2:1 OLR determination strategy presented in Trick et al. 
(2021). 

Wang et al. (2020) classified R — Va ridges into two 
types depending on whether they exhibit significant 
variation with stellar age. Coincidentally, the ridge dis- 
playing the most prominent variation with age roughly 
follows the constant Lz line at 2180kms~!kpc where 
Alz exhibits a peak feature. Given the tendency of 
older stars to be dynamically hotter, this prominent 
variation could be interpreted through the 2:1 outer 
Lindblad resonance of the Galactic Bar, from the dy- 
namic perspective. 


6. CONCLUSION 


Decoding the dynamical evolution history of the Milky 
Way is challenging due to the degeneracy caused by vari- 
ous perturbations. Our work attempts to shed new light 
on this task by differentiating the formation mechanism 
of the recently found phase space substructures within 
the Galactic Disk. The (Ve) corrugation pattern with 
Lz is the one-dimensional deprojection of R — Vg ridges 
and an easier target for comparative analysis. Leverag- 
ing the Gaia DR3 data and controlled numerical exper- 
iments using test particle simulations, we find that the 
Lz variation trend of phase shift between the Lz — (Vr) 
waves of different dynamical hotness could be an indi- 
cator of the dominant perturbations. With a simple toy 
model, our work also helps better comprehend the role of 
internal and external perturbers in shaping phase space 
substructures within the Galactic disk. 

Our key findings are summarized as follows: 

1. Analysis of the Gaia DR3 data reveals the Lz— (Vp) 
wave systematically shifts towards lower Lz in the dy- 
namically hotter population. Defining this phase shift 
as the difference in extrema location between the wave 
pattern of the cold and hot populations, its ampli- 
tude ALz exhibits a prominent peak feature at Lz ~ 
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2300 kms~!kpe which encodes information on the per- 
turbing mechanism. 

2. The external satellite perturbation produces a de- 
creasing trend with increasing Lz for ALz in the lower 
Lz range, which validates the theoretical prediction 
given by our toy model of radial phase mixing. This 
trend arises from the accumulation of radial frequency 
Opr difference between sub-populations of different dy- 
namical hotness. The persistence of this trend regard- 
less of the pericenter passages the satellite experienced 
makes it a unique signature of external perturbation. 

3. The corotating transient spirals do not produce a 
universal AL z variation trend with Lz for all azimuthal 
ranges. Combined with a steadily rotating bar, a char- 
acteristic ALz peak emerges at the 2:1 Outer Lindblad 
Resonance in most azimuthal ranges, which qualitatively 
resembles the observational feature. 

4. Comparisons between observation data and simu- 
lation results suggest an internal cause for the observed 
Lz — (Vr) wave, i.e. likely from bar and spiral arms. 
If that is the case, the Alz peak marks the location 
of 2:1 outer Lindblad resonance, which points towards a 
long/slow bar model with pattern speed between 33 and 
41 kms~tkpcu!. 

In conclusion, we demonstrate that examining the 
phase shift behaviors between dynamically hot and cold 
stellar populations holds the potential to discriminate 
between different types of perturbations. Additionally, 
by pinpointing a likely internal origin for the phase shift 
variation trends, our work provides constraints on the 
relative importance of internal and external processes in 
the recent dynamical evolution history of the Galactic 
disk. 

While our study represents significant progress, lim- 
itations exist that motivate several promising avenues 
for future investigations. For instance, tailored N-body 
simulations are required to investigate the self-gravity 


effects absent in test particle models. Reliable mathe- 
matical techniques must be developed for robustly mea- 
suring phase shifts, especially for non-sinusoidal wave 
shapes. Upcoming work can also examine phase shift 
variation patterns when combining different perturba- 
tions in test particle simulations. 

Future spectroscopic surveys such as 4MOST and 
Milky Way Mapper covering a wider range in spatial 
extent will be capable of uncovering the azimuthal vari- 
ation of phase shift variation pattern, which could help 
constrain possible perturbation scenarios. 


7. ACKNOWLEDGMENTS 


This work has made use of data from the Euro- 
pean Space Agency(ESA) mission Gaia (https://www. 
cosmos.esa.int/gaia), processed by the Gaia Data Pro- 
cessing and Analysis Consortium (DPAC, https://www. 
cosmos.esa.int/web/gaia/dpac/consortium). Funding 
for the DPAC has been provided by national institu- 
tions, in particular the institutions participating in the 
Gaia Multilateral Agreement. 

This work is supported by the National Natural Sci- 
ence Foundation of China under grant No.12122301, 
12233001, by a Shanghai Natural Science Research 
Grant (21ZR1430600), by the “111” project of the Min- 
istry of Education under grant No. B20019, and by the 
China Manned Space Project with No. CMS-CSST- 
2021-B03. We thank the sponsorship from Yangyang 
Development Fund. This work made use of the Grav- 
ity Supercomputer at the Department of Astronomy, 
Shanghai Jiao Tong University. 


Software: astropy (Astropy Collaboration et al. 
2013, 2018), agama(Vasiliev 2019), matplotlib(Hunter 
2007), galpy(Bovy 2015) 


Facilities: Gaia 


REFERENCES 


Anders, F., Khalatyan, A., Queiroz, A. B. A., et al. 2022, 
A&A, 658, A91, doi: 10.1051/0004-6361/202142369 

Antoja, T., Ramos, P., Lépez-Guitart, F., et al. 2022, 
A&A, 668, A61, doi: 10.1051/0004-6361/202244064 

Antoja, T., Helmi, A., Romero-Gémez, M., et al. 2018, 
Nature, 561, 360, doi: 10.1038/s41586-018-0510-7 

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., 
et al. 2013, A&A, 558, A33, 
doi: 10.1051/0004-6361/201322068 

Astropy Collaboration, Price-Whelan, A. M., Sipécz, B. M., 
et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881 /aabc4f 


Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417, 
doi: 10.1093/mnras/sty2813 

Binney, J. 2010, MNRAS, 401, 2318, 
doi: 10.1111/j.1365-2966.2009.15845.x 

—. 2012, MNRAS, 426, 1324, 
doi: 10.1111/j.1365-2966.2012.21757.x 

Binney, J., & McMillan, P. 2011, MNRAS, 413, 1889, 
doi: 10.1111/j.1365-2966.2011.18268.x 

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: 
Second Edition (Princeton University Press) 

Bovy, J. 2015, ApJS, 216, 29, 
doi: 10.1088 /0067-0049/216/2/29 


Carlberg, R. G., & Sellwood, J. A. 1985, ApJ, 292, 79, 
doi: 10.1086/163134 

Chakrabarty, D. 2007, A&A, 467, 145, 
doi: 10.1051/0004-6361:20066677 

Chiba, R., Friske, J. K. S., & Schénrich, R. 2021, MNRAS, 
500, 4710, doi: 10.1093/mnras/staa3585 

Chiba, R., & Schénrich, R. 2021, MNRAS, 505, 2412, 
doi: 10.1093/mnras/stab1094 

Cox, D. P., & Gómez, G. C. 2002, ApJS, 142, 261, 
doi: 10.1086/341946 

de la Vega, A., Quillen, A. C., Carlin, J. L., Chakrabarti, 
S., & D’Onghia, E. 2015, MNRAS, 454, 933, 
doi: 10.1093/mnras/stv2055 

Dehnen, W. 2000, AJ, 119, 800, doi: 10.1086/301226 

Fragkoudi, F., Katz, D., Trick, W., et al. 2019, MNRAS, 
488, 3324, doi: 10.1093/mnras/stz1875 

Friske, J. K. S., & Schénrich, R. 2019, MNRAS, 490, 5414, 
doi: 10.1093/mnras/stz2951 

Gaia Collaboration, Katz, D., Antoja, T., et al. 2018, A&A, 
616, A11, doi: 10.1051/0004-6361/201832865 

Gómez, F. A., Minchev, I., Villalobos, Á., O’Shea, B. W., 
& Williams, M. E. K. 2012, MNRAS, 419, 2163, 
doi: 10.1111/j.1365-2966.2011.19867.x 

GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 
2021, A&A, 647, A59, doi: 10.1051/0004-6361/202040208 

Hunt, J. A. S., Bub, M. W., Bovy, J., et al. 2019, MNRAS, 
490, 1026, doi: 10.1093/mnras/stz2667 

Hunt, J. A. S., Hong, J., Bovy, J., Kawata, D., & Grand, 
R. J. J. 2018, MNRAS, 481, 3794, 
doi: 10.1093/mnras/sty2532 

Hunter, J. D. 2007, Computing in Science & Engineering, 9, 
90, doi: 10.1109/MCSE.2007.55 

Kawata, D., Baba, J., Ciucă, I., et al. 2018, MNRAS, 479, 
L108, doi: 10.1093/mnrasl/sly107 


13 


Khanna, S., Sharma, S., Tepper-Garcia, T., et al. 2019, 
MNRAS, 489, 4962, doi: 10.1093/mnras/stz2462 
Laporte, C. F. P., Famaey, B., Monari, G., et al. 2020, 
A&A, 643, L3, doi: 10.1051/0004-6361/202038740 
Laporte, C. F. P., Minchev, I., Johnston, K. V., & Gómez, 
F. A. 2019, MNRAS, 485, 3134, 
doi: 10.1093/mnras/stz583 
Li, Z.-Y., & Shen, J. 2020, ApJ, 890, 85, 
doi: 10.3847 /1538-4357 /ab6b21 
Minchey, I., Quillen, A. C., Williams, M., et al. 2009, 
MNRAS, 396, L56, doi: 10.1111/j.1745-3933.2009.00661.x 
Monari, G., Famaey, B., Siebert, A., et al. 2016, MNRAS, 
461, 3835, doi: 10.1093/mnras/stw1564 
Monari, G., Famaey, B., Siebert, A., Wegg, C., & Gerhard, 
O. 2019, A&A, 626, A41, 
doi: 10.1051/0004-6361/201834820 
Miihlbauer, G., & Dehnen, W. 2003, A&A, 401, 975, 
doi: 10.1051/0004-6361:20030186 
Quillen, A. C., Carrillo, I., Anders, F., et al. 2018, MNRAS, 
480, 3132, doi: 10.1093/mnras/sty2077 
Ramos, P., Antoja, T., & Figueras, F. 2018, A&A, 619, 
A72, doi: 10.1051/0004-6361/201833494 
Reid, M. J., & Brunthaler, A. 2020, ApJ, 892, 39, 
doi: 10.3847 /1538-4357/ab76cd 
Trick, W. H. 2022, MNRAS, 509, 844, 
doi: 10.1093/mnras/stab2866 
Trick, W. H., Fragkoudi, F., Hunt, J. A. S., Mackereth, 
J. T., & White, S. D. M. 2021, MNRAS, 500, 2645, 
doi: 10.1093/mnras/staa3317 
Vasiliev, E. 2019, MNRAS, 482, 1525, 
doi: 10.1093/mnras/sty2672 
Wang, H. F., Huang, Y., Zhang, H. W., et al. 2020, ApJ, 
902, 70, doi: 10.3847/1538-4357 /abb3c8 


